cd "C:\Users\m1049\Dropbox\CPS Replication Final"

** It requires Stata 16 to produce Figures 3 and B4. 

use main.dta,replace

* Figure 3

gen Dyr=2010 if ct==6 | ct==10 |  ct==11 | ct==18
replace Dyr=2014 if ct==8

* watersupply

xtreg watersupply treatment_binary i.year, fe vce(cluster id)

preserve

gen D_4 = year - Dyr == -4
gen D_3 = year - Dyr == -3
gen D_2 = year - Dyr == -2
gen D_1 = year - Dyr == -1
gen D0 = year - Dyr == 0
gen D1 = year - Dyr == 1
gen D2 = year - Dyr == 2
gen D3 = year - Dyr == 3
gen D4 = year - Dyr == 4
gen D5 = year - Dyr == 5
gen D6 = year - Dyr == 6
gen D7 = year - Dyr == 7
gen D8 = year - Dyr == 8
gen D9 = year - Dyr == 9

xtreg watersupply D_4 D_3 D_2 D0 D1 D2 D3 D4 D5 D6 D7 D8 D9 i.year, fe vce(cluster id)

parmest, fast level(95 90)
gen parmseq=_n

keep if parmseq <= 13
replace parmseq = parmseq - 4 if parmseq>=4
replace parmseq = parmseq - 5 in 1/3
set obs 14
replace parm = "D_1" in 14
replace parmseq = -1 in 14
replace estimate = 0 in 14

#delimit;

gr tw

	(rcap min95 max95 parmseq, lpattern(solid) lcolor(gs8))
	(scatter est parmseq, msymbol(o) msize(vlarge) mcolor(black))

	,
		yline(0, lpat(solid) lwid(vthin))
		xline(0, lpat(dash) lcolor(blue))
		xlab(-4(1)9)
		plotregion(style(none))
		legend(off)
		ylab(, angle(horizontal) format(%10.0f))
		l2title("Difference (Treated-Control)")
		scheme(s1mono)
		xtitle("Years since the switch")
		ytitle("") 
		title("Tap Water Coverage")
		saving(leads_water,replace)
		;
#delimit cr

restore

* road

xtreg road treatment_binary i.year, fe vce(cluster id)

preserve

gen D_4 = year - Dyr == -4
gen D_3 = year - Dyr == -3
gen D_2 = year - Dyr == -2
gen D_1 = year - Dyr == -1
gen D0 = year - Dyr == 0
gen D1 = year - Dyr == 1
gen D2 = year - Dyr == 2
gen D3 = year - Dyr == 3
gen D4 = year - Dyr == 4

xtreg road D_4 D_3 D_2 D0 D1 D2 D3 D4 i.year, fe vce(cluster id) 

parmest, fast level(95 90)

keep in 1/8
gen parmseq=.
replace parmseq =- 4 in 1
replace parmseq =- 3 in 2
replace parmseq =- 2 in 3
replace parmseq = 0 in 4
replace parmseq = 1 in 5
replace parmseq = 2 in 6
replace parmseq = 3 in 7
replace parmseq = 4 in 8

set obs 9
replace parm = "D_1" in 9
replace parmseq = -1 in 9
replace estimate = 0 in 9

#delimit;

gr tw

	(rcap  min95 max95 parmseq, lpattern(solid) lcolor(gs8))
	(scatter est parmseq, msymbol(o) msize(vlarge) mcolor(black))

	,
		yline(0, lpat(solid) lwid(vthin))
		xline(0, lpat(dash) lcolor(blue))
		xlab(-4(1)4)
		plotregion(style(none))
		legend(off)
		ylab(, angle(horizontal) format(%10.0f))
		l2title("Difference (Treated-Control)")
		xtitle("Years since the switch")
		title("Motor Vehicle Fatality Rate")
		scheme(s1mono)
		ytitle("") saving(leads_road,replace)
		;
#delimit cr

restore

* library

xtreg lib treatment_binary i.year, fe vce(cluster id)

preserve 

gen D_2 = year - Dyr == -2
gen D_1 = year - Dyr == -1
gen D0 = year - Dyr == 0
gen D1 = year - Dyr == 1
gen D2 = year - Dyr == 2
gen D3 = year - Dyr == 3
gen D4 = year - Dyr == 4
gen D5 = year - Dyr >= 5

xtreg lib D_2 D0 D1 D2 D3 D4 i.year, fe vce(cluster id)

parmest, fast level(95 90)

keep in 1/6
gen parmseq=.
replace parmseq =-2 in 1
replace parmseq =0 in 2
replace parmseq =1 in 3
replace parmseq =2 in 4
replace parmseq =3 in 5
replace parmseq =4 in 6

set obs 7
replace parm = "D_1" in 7
replace parmseq = -1 in 7
replace estimate = 0 in 7

#delimit;

gr tw

	(rcap min95 max95 parmseq, lpattern(solid) lcolor(gs8))
	(scatter est parmseq, msymbol(o) msize(vlarge) mcolor(black))

	,
		yline(0, lpat(solid) lwid(vthin))
		xline(0, lpat(dash) lcolor(blue))
		xlab(-2(1)4)
		plotregion(style(none))
		legend(off)
		ylab(, angle(horizontal) format(%10.0f))
		l2title("Difference (Treated-Control)")
		xtitle("Years since the switch")
		scheme(s1mono)
		ytitle("") 
		title("Community Libraries")
		saving(leads_lib,replace)
		;
#delimit cr

restore 

* activity center

xtreg center treatment_binary i.year, fe vce(cluster id)

preserve 

gen D_2 = year - Dyr == -2
gen D_1 = year - Dyr == -1
gen D0 = year - Dyr == 0
gen D1 = year - Dyr == 1
gen D2 = year - Dyr == 2
gen D3 = year - Dyr == 3
gen D4 = year - Dyr == 4

xtreg center D_2 D0 D1 D2 D3 D4 i.year, fe vce(cluster id)

parmest, fast level(95 90)

keep in 1/6
gen parmseq=.
replace parmseq =-2 in 1
replace parmseq =0 in 2
replace parmseq =1 in 3
replace parmseq =2 in 4
replace parmseq =3 in 5
replace parmseq =4 in 6

set obs 7
replace parm = "D_1" in 7
replace parmseq = -1 in 7
replace estimate = 0 in 7

#delimit;

gr tw

	(rcap min95 max95 parmseq, lpattern(solid) lcolor(gs8))
	(scatter est parmseq, msymbol(o) msize(vlarge) mcolor(black))
	,
		yline(0, lpat(solid) lwid(vthin))
		xline(0, lpat(dash) lcolor(blue))
		xlab(-2(1)4)
		plotregion(style(none))
		legend(off)
		ylab(, angle(horizontal) format(%10.0f))
		l2title("Difference in Community Centers")
		xtitle("Years since the switch", height(4))
		scheme(s1mono)
		ytitle("") 
		title("Community Centers")
		saving(leads_center,replace)
		;
#delimit cr

restore

* storage

xtreg max treatment_binary i.year, fe vce(cluster id)

preserve

gen D_2 = year - Dyr == -2
gen D_1 = year - Dyr == -1
gen D0 = year - Dyr == 0
gen D1 = year - Dyr == 1
gen D2 = year - Dyr == 2
gen D3 = year - Dyr == 3

xtreg max D_2 D0 D1 D2 D3 i.year , fe vce(cluster id)

parmest, fast level(95 90)

keep in 1/5
gen parmseq=.
replace parmseq =-2 in 1
replace parmseq =0 in 2
replace parmseq =1 in 3
replace parmseq =2 in 4
replace parmseq =3 in 5

set obs 6
replace parm = "D_1" in 6
replace parmseq = -1 in 6
replace estimate = 0 in 6

#delimit;

gr tw

	(rcap min95 max95 parmseq, lpattern(solid) lcolor(gs8))
	(scatter est parmseq, msymbol(o) msize(vlarge) mcolor(black))

	,
		yline(0, lpat(solid) lwid(vthin))
		xline(0, lpat(dash) lcolor(blue))
		xlab(-2(1)3)
		plotregion(style(none))
		legend(off)
		ylab(, angle(horizontal) format(%10.0f))
		l2title("Difference (Treated-Control)")
		xtitle("Years since the switch")
		scheme(s1mono)
		ytitle("") 
		title("Max Storage")
		saving(leads_storage,replace)
		;
#delimit cr

restore 

graph combine leads_water.gph leads_road.gph leads_lib.gph leads_center.gph leads_storage.gph, scheme(s1mono) note("95% Confidence Intervals" "The omitted period is t = −1")

graph export leads_combine.pdf,replace

* Figure 4

use main.dta,replace

* water

xtreg water treatment_binary##c.elite_capture i.year, fe vce(cluster id)

frame create mplot vote_buying coefficient lower upper

forvalues a = -0.25(0.15)0.9 {

lincom _b[1.treatment_binary] + _b[1.treatment_binary#c.elite_capture]*`a'

frame post mplot (`a') (`r(estimate)')(`r(lb)') (`r(ub)')

}

frame change mplot
isid vote_buying, sort
format coefficient %9.0g
format lower %9.0g
format upper %9.0g
list, noobs clean

graph twoway  (rspike lower upper vote_buying) (scatter coefficient vote_buying , sort), scheme(s1mono) yline(0,lpattern(dash)) xscale(range (-0.25 0.85)) xtitle("Latent Measure of Local Capture", height(5)) legend(off) ytitle("Predicted Conditional DiD Estimate", height(5)) title("Tap Water", margin(b=2.5)) saving(mplot_water_capture,replace)

frame change default
frame drop mplot

* road

xtreg road treatment_binary##c.elite_capture i.year, fe vce(cluster id)

frame create mplot vote_buying coefficient lower upper

forvalues a = -0.25(0.15)0.9 {

lincom _b[1.treatment_binary] + _b[1.treatment_binary#c.elite_capture]*`a'

frame post mplot (`a') (`r(estimate)')(`r(lb)') (`r(ub)')

}

frame change mplot
isid vote_buying, sort
format coefficient %9.0g
format lower %9.0g
format upper %9.0g
list, noobs clean

graph twoway  (rspike lower upper vote_buying) (scatter coefficient vote_buying , sort), scheme(s1mono) yline(0,lpattern(dash)) xscale(range (-0.25 0.85)) xtitle("Latent Measure of Local Capture", height(5)) legend(off) ytitle("Predicted Conditional DiD Estimate", height(5)) title("Road Conditions", margin(b=2.5)) saving(mplot_road_capture,replace)

frame change default
frame drop mplot

* lib

xtreg lib treatment_binary##c.elite_capture i.year, fe vce(cluster id)

frame create mplot vote_buying coefficient lower upper

forvalues a = -0.25(0.15)0.9 {

lincom _b[1.treatment_binary] + _b[1.treatment_binary#c.elite_capture]*`a'

frame post mplot (`a') (`r(estimate)')(`r(lb)') (`r(ub)')

}

frame change mplot
isid vote_buying, sort
format coefficient %9.0g
format lower %9.0g
format upper %9.0g
list, noobs clean

graph twoway  (rspike lower upper vote_buying) (scatter coefficient vote_buying , sort), scheme(s1mono) yline(0,lpattern(dash)) xscale(range (-0.25 0.85)) xtitle("Latent Measure of Local Capture", height(5)) legend(off) ytitle("Predicted Conditional DiD Estimate", height(5)) title("Community Libraries", margin(b=2.5)) saving(mplot_lib_capture,replace)

frame change default
frame drop mplot

* center

xtreg center treatment_binary##c.elite_capture i.year, fe vce(cluster id)

frame create mplot vote_buying coefficient lower upper

forvalues a = -0.25(0.15)0.9 {

lincom _b[1.treatment_binary] + _b[1.treatment_binary#c.elite_capture]*`a'

frame post mplot (`a') (`r(estimate)')(`r(lb)') (`r(ub)')

}

frame change mplot
isid vote_buying, sort
format coefficient %9.0g
format lower %9.0g
format upper %9.0g
list, noobs clean

graph twoway  (rspike lower upper vote_buying) (scatter coefficient vote_buying , sort), scheme(s1mono) yline(0,lpattern(dash)) xscale(range (-0.25 0.85)) xtitle("Latent Measure of Local Capture", height(5)) legend(off) ytitle("Predicted Conditional DiD Estimate", height(5)) title("Community Centers", margin(b=2.5)) saving(mplot_center_capture,replace)

frame change default
frame drop mplot

* max storage

xtreg max treatment_binary##c.elite_capture i.year, fe vce(cluster id)

frame create mplot vote_buying coefficient lower upper

forvalues a = -0.25(0.15)0.9 {

lincom _b[1.treatment_binary] + _b[1.treatment_binary#c.elite_capture]*`a'

frame post mplot (`a') (`r(estimate)')(`r(lb)') (`r(ub)')

}

frame change mplot
isid vote_buying, sort
format coefficient %9.0g
format lower %9.0g
format upper %9.0g
list, noobs clean

graph twoway (rspike lower upper vote_buying) (scatter coefficient vote_buying, sort), scheme(s1mono) yline(0,lpattern(dash)) xscale(range (-0.25 0.85)) xtitle("Latent Measure of Local Capture", height(5)) legend(off) ytitle("Predicted Conditional DiD Estimate", height(5)) title("Max Storage", margin(b=2.5)) saving(mplot_max_capture,replace)

frame change default
frame drop mplot

graph combine mplot_road_capture.gph mplot_lib_capture.gph mplot_center_capture.gph mplot_max_capture.gph, scheme(s1mono) 

graph export hte_capture_combine.pdf,replace

* Figure B.1

foreach var of varlist  mean_income pan_blue_legis_vote population_density net_migration primary_education elite_capture{
	
bysort id: egen premean_`var'=mean(`var') if year<=2009 & year>=2005

}


foreach var of varlist premean_mean_income premean_pan_blue_legis_vote premean_population_density premean_net_migration premean_primary_education premean_elite_capture{

egen mean`var'=mean(`var')
egen sd`var' = sd(`var')
gen z`var'=(`var'-mean`var')/sd`var'
}

reg zpremean_elite_capture zpremean_pan_blue_legis_vote zpremean_primary_education zpremean_population_density zpremean_mean_income zpremean_net_migration if year==2009, r

label variable zpremean_mean_income "Mean Income"
label variable zpremean_pan_blue_legis_vote "KMT Vote Share"
label variable zpremean_population_density "Population Density"
label variable zpremean_net_migration "Net Migration Rate"
label variable zpremean_primary_education "Primary Education"

coefplot, drop(_cons) xline(0, lpattern(dash) lcolor(gray)) levels(90) scheme(s1mono)

* Figure B.3

* road

use "main", replace

xtreg road treatment_binary i.year, fe vce(cluster id)

keep if e(sample)

egen id1=group(id)

global tflist ""
                forval i = 1/303 {
					tempfile tfcur
					xtreg road treatment_binary i.year if id1!=`i', fe vce(cluster id)

					parmest, label format(estimate min95 max95 %8.2f p %8.1e) idn(`i') saving(`"`tfcur'"',replace) flist(tflist)
				}
				
drop _all		
append using $tflist
sort idnum 

keep  idnum parm label  estimate p
encode parm, gen (variable)
label list variable

keep if variable==17

sort idnum
drop parm label
reshape wide estimate p, i(idnum) j(variable)

rename estimate17 DiD
rename p17 pvalue_DiD

scatter pvalue_DiD DiD, xtitle("DiD", height(5)) ytitle("P-Value", height(5)) ylabel(0 0.00004 0.00008 ,format(%9.0g)) xlabel( ,format(%9.0g)) scheme(s1mono) title("Roads", margin(b=2.5)) saving("roads",replace)

* lib

use "main.dta",replace

xtreg lib treatment_binary i.year, fe vce(cluster id)
keep if e(sample)

egen id1=group(id)


global tflist ""
                forval i = 1/133 {
					tempfile tfcur
					xtreg lib treatment_binary i.year if id1!=`i', fe vce(cluster id)

					parmest, label format(estimate min95 max95 %8.2f p %8.1e) idn(`i') saving(`"`tfcur'"',replace) flist(tflist)
				}
				
drop _all		
append using $tflist
sort idnum 

keep  idnum parm label  estimate p
encode parm, gen (variable)
label list variable

keep if variable==9

sort idnum
drop parm label
reshape wide estimate p, i(idnum) j(variable)

rename estimate9 DiD
rename p9 pvalue_DiD

drop if idnum==124 | idnum==106
gsort -pvalue_DiD

scatter pvalue_DiD DiD, xtitle("DiD", height(5)) ytitle("P-Value", height(5)) ylabel( ,format(%9.0g)) xlabel( ,format(%9.0g)) scheme(s1mono) title("Libraries", margin(b=2.5)) saving("lib",replace)

* center

use "main.dta",replace

xtreg center treatment_binary i.year, fe vce(cluster id)
keep if e(sample)

egen id1=group(id)

sum id1

global tflist ""
                forval i = 1/139 {
					tempfile tfcur
					xtreg center treatment_binary i.year if id1!=`i', fe vce(cluster id)

					parmest, label format(estimate min95 max95 %8.2f p %8.1e) idn(`i') saving(`"`tfcur'"',replace) flist(tflist)
				}
				
drop _all		
append using $tflist
sort idnum 

keep  idnum parm label  estimate p
encode parm, gen (variable)
label list variable

keep if variable==9

sort idnum
drop parm label
reshape wide estimate p, i(idnum) j(variable)

rename estimate9 DiD
rename p9 pvalue_DiD

scatter pvalue_DiD DiD, xtitle("DiD", height(5)) ytitle("P-Value", height(5)) ylabel( ,format(%9.0g)) xlabel( ,format(%9.0g)) scheme(s1mono) title("Centers", margin(b=2.5)) saving("centers",replace)

* gray

use "main.dta",replace

xtreg max treatment_binary i.year, fe vce(cluster id)

keep if e(sample)

egen id1=group(id)

global tflist ""
                forval i = 1/226 {
					tempfile tfcur
					xtreg max treatment_binary i.year if id1!=`i', fe vce(cluster id)

					parmest, label format(estimate min95 max95 %8.2f p %8.1e) idn(`i') saving(`"`tfcur'"',replace) flist(tflist)
				}
				
drop _all		
append using $tflist
sort idnum 

keep  idnum parm label  estimate p
encode parm, gen (variable)
label list variable

keep if variable==8

sort idnum
drop parm label
reshape wide estimate p, i(idnum) j(variable)

rename estimate8 DiD
rename p8 pvalue_DiD

scatter pvalue_DiD DiD, xtitle("DiD", height(5)) ytitle("P-Value", height(5)) ylabel( ,format(%9.0g)) xlabel( ,format(%9.0g)) scheme(s1mono) title("Columbarium Niches", margin(b=2.5)) saving("niches",replace)

graph combine roads.gph lib.gph centers.gph niches.gph, row(2) scheme(s1mono)
graph export jackknife_combine.pdf,replace

* Figure B.4

use main.dta,replace

* water

xtreg water treatment_binary##c.primary i.year, fe vce(cluster id)
frame create mplot primary coefficient lower upper
forvalues a = 0.2(0.1)0.8 {
    lincom _b[1.treatment_binary] + _b[1.treatment_binary#c.primary]*`a'
    frame post mplot (`a') (`r(estimate)')(`r(lb)') (`r(ub)')
}

frame change mplot
isid primary, sort
format coefficient %9.0g
format lower %9.0g
format upper %9.0g
list, noobs clean

graph twoway  (rspike lower upper primary) (scatter coefficient primary , sort), scheme(s1mono) yline(0,lpattern(dash)) xscale(range (0.15 0.85)) xtitle("Primary Education below (share)", height(5)) legend(off) ytitle("DiD Estimate", height(5)) title(Tap Water Supply,margin(b=2.5)) saving(mplot_water_primary,replace)

frame change default
frame drop mplot

* road

xtreg road treatment_binary##c.primary i.year, fe vce(cluster id)
frame create mplot primary coefficient lower upper
forvalues a = 0.2(0.1)0.8 {
    lincom _b[1.treatment_binary] + _b[1.treatment_binary#c.primary]*`a'
    frame post mplot (`a') (`r(estimate)')(`r(lb)') (`r(ub)')
}

frame change mplot
isid primary, sort
format coefficient %9.0g
format lower %9.0g
format upper %9.0g
list, noobs clean

graph twoway  (rspike lower upper primary) (scatter coefficient primary , sort), scheme(s1mono) yline(0,lpattern(dash)) xscale(range (0.15 0.85)) xtitle("Primary Education below (share)", height(5)) legend(off) ytitle("DiD Estimate", height(5)) title(Road Conditions, margin(b=2.5)) saving(mplot_road_primary,replace)

frame change default
frame drop mplot

* lib

xtreg lib treatment_binary##c.primary i.year, fe vce(cluster id)
frame create mplot primary coefficient lower upper
forvalues a = 0.2(0.1)0.8 {
    lincom _b[1.treatment_binary] + _b[1.treatment_binary#c.primary]*`a'
    frame post mplot (`a') (`r(estimate)')(`r(lb)') (`r(ub)')
}

frame change mplot
isid primary, sort
format coefficient %9.0g
format lower %9.0g
format upper %9.0g
list, noobs clean

graph twoway  (rspike lower upper primary) (scatter coefficient primary , sort), scheme(s1mono) yline(0,lpattern(dash)) xscale(range (0.15 0.85)) xtitle("Primary Education below (share)", height(5)) legend(off) ytitle("DiD Estimate", height(5)) title(Community Libraries, margin(b=2.5)) saving(mplot_lib_primary,replace)

frame change default
frame drop mplot

* center

xtreg center treatment_binary##c.primary i.year, fe vce(cluster id)
frame create mplot primary coefficient lower upper
forvalues a = 0.2(0.1)0.8 {
    lincom _b[1.treatment_binary] + _b[1.treatment_binary#c.primary]*`a'
    frame post mplot (`a') (`r(estimate)')(`r(lb)') (`r(ub)')
}

frame change mplot
isid primary, sort
format coefficient %9.0g
format lower %9.0g
format upper %9.0g
list, noobs clean

graph twoway  (rspike lower upper primary) (scatter coefficient primary , sort), scheme(s1mono) yline(0,lpattern(dash)) xscale(range (0.15 0.85)) xtitle("Primary Education below (share)", height(5)) legend(off) ytitle("DiD Estimate", height(5)) title(Community Centers, margin(b=2.5)) saving(mplot_center_primary,replace)

frame change default
frame drop mplot

* max

xtreg max treatment_binary##c.primary i.year, fe vce(cluster id)
frame create mplot primary coefficient lower upper
forvalues a = 0.2(0.1)0.8 {
    lincom _b[1.treatment_binary] + _b[1.treatment_binary#c.primary]*`a'
    frame post mplot (`a') (`r(estimate)')(`r(lb)') (`r(ub)')
}

frame change mplot
isid primary, sort
format coefficient %9.0g
format lower %9.0g
format upper %9.0g
list, noobs clean

graph twoway  (rspike lower upper primary) (scatter coefficient primary , sort), scheme(s1mono) yline(0,lpattern(dash)) xscale(range (0.15 0.85)) xtitle("Primary Education below (share)", height(5)) legend(off) ytitle("DiD Estimate", height(5)) title(Max Storage, margin(b=2.5)) saving(mplot_max_primary,replace)

frame change default
frame drop mplot

graph combine mplot_water_primary.gph mplot_road_primary.gph mplot_lib_primary.gph mplot_center_primary.gph mplot_max_primary.gph, scheme(s1mono) 

graph export hte_primary_combine.pdf,replace
